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PREFACE 


This Annual Technical Report presents research performed during the 
first year of NASA-Lewls Grant NSC-3217 • which was Initiated on August 1, 
1978. The grant is continuing, presently being budgeted for the second 
year. The NASA-Lewls Technical Monitor since the Inception of this grant 
has been Dr. J. A. DlCarlo the Materials Science Branch. 

The study Is being performed within the Composite Materials Research 
Group at the University of Wyoming. Principal Investigators are Mr. 

Daniel P. Murphy, Graduate Student, and Dr. Donald F. Adams, Professor, 
both of the Mechanical Engineering Department. Mohamed M. Monlb and Brent 
G. Schaffer, Graduate Students In Mechanical Engineering, have also made 
significant contributions. 
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SECTION I 


INTRODUCTION 


It is widely recognized that crack initiation and propagation in 
fiber-reinforced composite materials is necessarily quite different from 
that in homogenous solids, and that most applications of classical lin- 
early elastic fracture mechanics have not been satisfactory. 

The Inhomogenous nature of fiber-reinforced materials strongly 
suggests the use of so-called micromechanical analysis techniques, where- 
in the interactions of the individual reinforcing fibers with the surround 

Ik. 

Ing matrix material are considered. Specifically, if sufficient insight 
into the geometrical effects can be gained, the energy absorption due to 
crack initiation and growth in a unidirectional boron/ aluminum composite 
can be characterized with a two-dimensional micromechanics analysis. 

This Is the primary objective of the present study. 

The first-year effort reported here has been directed at investigat- 
ing the response of a unidirectional boron/aluminum composite containing 
a fiber flaw site or discontinuity. Tensile loading in the direction of 
the fiber axes is of primary Interest, although transverse loading has 
also been analyzed. To approximate the behavior of a metal matrix com- 
posite under these conditions, two features are of major importance to a 
micromechanics analysis : 

“The ability to model the full elastic-plastic range of the 
matrix material and modify it to account for changes in 
temperature, 

“The availability of a procedure for approximating crack 


growth in the matrix materiol. 

The finite element micromechanica program previoualy developed at 
the University of Wyoming (1] possessed the first of these features at 
the time that the present grant was awarded, and a crack propagation 
capability was incorporated during the fourth quarter of this first- 
year effort. 

A survey of the literature available on the topics of mlcromechan- 
Ics analyses and the effect of flaws in unidirectional metal matrix com- 
posites has yielded some valuable information, and has provided insight 
into the problems involved. But the exact problem dealt with in this 
grant study does not appear to have been attempted before. 

Among previous analysis efforts is a series of two reports published 
in 1973 by Adams [2], and Repnau and Adams [3], in which a finite ele- 
ment micromechanics program was developed, incorporating both elasto- 
plastic capability and crack propagation. A unidirectional boron/aluminum 
composite was studied, although no flaw sites were incorporated and the 
model grid represented a section of the composite which was transverse 
to the fiber axes. Also, because a plane strain solution was assumed, 
only transverse loading could be studied. Material modifications to 
reflect environmental changes were also not available in these programs. 

Another effort was that of Ko [4] In 1977. His analysis was of an 
axially loaded boron/aluminum composite using the NASTRAN finite element 
analysis program as a micromechanics tool. His model consisted of a 
single boron fiber surrounded by an annular section of aluminum matrix. 
Axlsymmetrlc finite elements were employed and discontinultltes In both 
the fiber and the matrix were studied. However, the NASTRAN solution did 
not Incorporate elastoplastlc material response, nor was any sort of 
crack propagation scheme used. 
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Akbarzadeh [5] undertook a mlcroroechanlca analysis of flawed uni- 
directional composites In 1976. The material he studied was E-glasa/ 
epoxy and the flaws were taken to be discontinuous fibers. Again, the 
NASTRAN finite element program was employed, and no elastoplastic de- 
formation or crack growth was considered. His finite element grid rep- 
resented a cross-sectional cut transverse to the fiber axes and all load- 
ing was confined to the plane that section, i.e. , there was no axial 
loading. Akbarzadeh's finite elemem; grid was very fine, and his analy- 
sts yielded valuable information concerning the elastic microstress state 
in composites with different fiber packing geometries and densities, and 
the effects of varying the clastic moduli of the constituent materials. 

Finally, there has been a limited amount of experimental work done 
in the area of interest dealt with in this report. One of the more im- 
portant of these was by Awerbuch and Hahn [6], in which the crack tip dam- 
age and fracture toughness of axially lo.aded unidirectional boron/aluminum 
composites was stvidled. Awerbuch and Hahn prepared tension coupons of 
boron/aluminum in which a center notch was machined, cutting several fibers 
and the aluminum matrix between them. Among the data generated by their 
efforts are some describing the deformation of the aluminum matrix surround- 
ing the last cut fiber at the edge of the notch. This has been of some 
help in evaluating the results of the present analysis. 

The present study has yielded valuable information about the analyti- 
cal approach to the problem of energy absorption in a flawed metal matrix 
composite. The analysis is a quasi-three-dlmensional formulation (see 
Section 3), which allows traction loads to be applied to threee mutually 
perpendicular surfaces ot a two-dimensional finite element array. This 
program in its present form can predict the microstress state of a metal 
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macrix compoaite wlch flawa quite veil, and providee infonxation regard- 
ing the energy abaorption reaulting from plaatic deformation and crack 
growth at auch flow aitea. However, the full characterization of thia 
problem involvea the rather complex geometry reaulting from on array of 
cylindrical fibera inbedded in an elaatoplaatic matrix, and the final 
aolutlon would appear to require a three-dimensional analysis. The follow- 
on second-year study will consider such an analysis, Incorporating all 
of the features of the present analysis into a three-dimensional form- 


ulation. 


SECTION 2 


SUMMARY 

Consldernble progress In the invcBtlgation of the energy absorption 
and the mlcrostress state of unidirectional metal matrix composites has 
been made to date. Further development of the analytical procedures de- 
scribed In this report, together with a carefully conceived and executed 
experimental program is expected to lead to a reliable means of pre- 
dicting the strength and toughness of these composite materiiils under 
actual conditions of manufacture and service. That is, given certain 
statistical parameters which define the quality of the constituents, a 
prediction of the energy capacity of the composite could be made. 

As was discussed in Section 1, and confirmed by an on-going litera- 
ture survey, the particular problem being dealt with here, and the analyti- 
cal procedures employed, make this study unique. A finite element analy- 
sis program is used to determine the stress-strain state of a unidirection- 
al, metal matrix composite material. A rather comprehensive description 
of the theoretical foundations and the special capabilities of this computer 
program Is provided in Appendix A. This program has been modified, and 
is in the process of being expanded to make it more suitable for the study 
of energy absorption mechanisms in axially loaded metal matrix composite 
materials containing flaws in Che form of microcracks in the matrix 
material, or discontinuities in the fibers. The two most Important energy 
absorption mechanisms being studied are plastic deformation due to the 
stress concentrations arising out of a material flaw, and growth of any 
cracks Initiated by such a flaw. 


A major problem in Invearlgntlng the mlcroBtrette state of uni- 
directional compoHltoM 1 b related to their geometrical Inhomogenelty, l.e., 
the fact that a detailed repreventatlon of the streBS state In and around 
the filaments embedded In an Inelastic matrix Is required. While the 
finite element analysis program employed Is quasl-three-dlmensional In 
that it can be loaded In all three coordinate directions and the result- 
ing stresses and strains computed, the actual modeling of the region of 
Interest Is essentially a two-dimensional representation. As described 
in Sections 3.3 and 3.4, an attempt to account for the complex geometry 
of the problem was made by constructing two types of finite element 
arrays. One Is the "longitudinal" model, which represents a cross-section- 
al cut through the composite running parallel to the ixes of the boron 
fibers. This model Is particularly useful in studying the effects of dis- 
continuous fibers, or fibers with regions of reduced strength. In addition, 
the presence of cracks or voids In the metal matrix can be modeled and their 
effect on the axial strength, and to some extent the transverse strength, 
of the material can be studied. 

The other type of model developed will be referred to the "transverse" 
section model, and depicts a cross-sectional cut perpendicular to the 
fiber axes of the composite. This model Is well-suited for studying the 
effects of transverse loading of the composite material, and the Influence 
of the circular cross section of the fibers is fully accounted for. Fiber 
flaws and matrix flaws can be completely characterized for transverse 
loading. Due to the generalized plane strain formulation used In the 
analysis. It Is possible to load the transverse model in the axial, or 
z-directlons. Unfortunately, with this mode of loading, it has not been 
possible to study to stress concentrating effects of flaws and dlscon- 
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tlnuleiaa. Thla la due to the fact that all element dlsplacementa out 
of the plane of the model muat be equal (ace Sect^oa A*2 of the Appendix), 
and variatlona in axial atreaa are thua due to Poiai?on ei^fecte and differ* 
eiicea in material propertiea onl}^ 

The longitudinal models have been useful in studying the effects of 
flaws on the axial strength of unidirectional composites, as detailed in 
Sections 5.1 and 5.2. The extent of plastic deformation around a flaw 
and the effects on the load redistribution among the surrounding fibers 
can be characterized quite well. For this model too, considerations 
have to be made regarding the geometry of the composite. As discussed 
in Section 3.3, this leads to the two configurations of the longitudinal 
model which are necessary to characterize a square or rectangluar fiber 
array. 

Crack initiation and propagation are clearly important considerations 
in the study of energy absorption in flawsd materials. This capability 
has been added to the micromechanics analysis, and is described in Section 
3.2. Unfortunately, a numerical inconsistency was discovered near the 
end of this first-year program, and subsequent corrections to the program 
have not yet been fully debugged. Consequently, numerical results of 
crack propagation in a unidirectional, metal matrix composite will be 
presented in the next quarterly report. 


SECTION 3 


ANALYSIS METHOD 

3.1 THE FINITE ELEMENT METHOD AS A MICROMECHANICS TOOL 

The scope of the following analyeis le to determine the internal 
stresB distribution in a unidirectional fiber-reinforced, metal matrix 
composite material subjected to mechanical loadings and variations in 
temperature. In addition, the effects of mater'lnl flaws or discontin- 
uities and their subsequent propagation through the material continuum 
are to be characterized. 

The typical unidirectional composite has a nonhomogenous internal 
structure that cdnsists of at least two distinct phases, i.e,, a homogen- 
ous matrix material reinforced by isotropic or transversely isotropic 
fibers. Transverse isotropy refers to the condition in which the axial 
properties of the fiber (e.g., strength and modulus) differ significantly 
from chose in a plane normal to the fiber axis, in most unidirectional 
composites, boron /aluminum Included, the reinforcing fibers are much 
stronger and possess a much higher axial elastic modulus chan the surround- 
ing matrix. The nonhomogenous nature of such a composite, together with 
the geometrical considerations of a cylindrical filament embedded in a 
matrix and containing a microscopic flaw, results in a boundary value 
problem of such compexity that a classical closed form ccntinuum solution 
to evaluate the microstress state would be impractical. As a result, in- 
vestigators for the last 12 to lA years have formulated numerical schemes 
to evaluate mlcrostresses in unidirectional composites. In one of the 
earlier efforts, Adams and Doner [7] applied finite difference techniques 
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to thttlr micronachonlcal atudioa. Kora racaiitly, the finite eleaent 
mthod haa amargod aa the moat varaatile tool available for aueh atudlea, 
aa denonatratad by nuonroua Invaatlgatora {8-18). 

The finite element analyala B»thod la baaed on the concept of dla- 
cretlzlng a material continuum Into an aaaembly of "elementa'* of "finite" 
alae. The Individual elementa are aaaemblad Into a network repreaentlng 
the continuum by joining them at predetermined polnta or "nodea" along 
their boundarlea. For any element • approximate functlona repreaentlng 
either atreaa, atraln or the dlaplacement field within that element can 
then be written. By a suitable choice of coefficients for these assumed 
field equations, a minimization of the potential energy of the system Is 
achieved. To dat(^ , the assumed displacement field technique has been the 
most successful and the most generally applied {191. For an assumed dis- 
placement field for the Interior of an assembly of elements, minimization 
of the potential energy results In a set of simultaneous algebraic equa- 
tions relating loads to displacements. These equations can readily be 
solved with modern digital computers. This procedure is expli.lned in 
rigorous detail by Heubner [19] and Zienklewlcz [20] v;herein they point 
out that the finite element method Is as useful in solving thermodynamics 
and fluid mechanics problems as It Is In solid mechanics, the only differ- 
ence being that functionals other than potential energy must be minimized. 

One of the great advantages of this method is that the approximate 
field equations need only satisfy the constraints of the individual ele- 
ments. An Important consequence of this situation is that equilibrium and 
compatibility conditions between the assumed fields of the individual ele- 
ments must be met in order to insure convergence to a correct solution. 
These conditions are also thoroughly discussed and developed in texts by 
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Heubner [19], Cook [21], and Zlcnkiewlcz [20]. 

The finite element method may be applied to any general three-di- 
mensional problem, the primary limitation being one of computer capacity. 

In an effort to analyze rather complex structural configurations In suffi- 
cient detail, most analyses. Including all micromechanlcs analyses to date, 
are reduced to two-dimensional problems, i.e., plane stress, plane strain, 
or generalized plane strain formulations. The analysis method used In the 
present study is a two-dimensional form Incorporating a condition of gen- 
eralized plane strain, which permits a specific loading to be applied in 
the third direction. This concept Is discussed In detail In Section A. 2 
of Appendix A. 

3.2 FINITE ELEMENT MICROMECHANICS ANALYSIS 

The primary analytical tool used in the present study has been the 
micromechanlcs finite element analysis program developed by A. K. Miller 
while completing his Ph.D. studies at the University of Wyoming, which Is 
described In considerable detail In Reference [1]. This program was created 
to Investigate the microstress state in unidirectional composite materials 
subjected to transverse mechanical loads, thermal gradients, and dllatatlon- 
al stresses due to moisture absorption by polymeric matrix materials. 

Among the special features of this program are Its ability to model the 
elastic-plastic stress-strain response of the isotropic matrix material, 
and in concert with the determination of thermal and moisture dilatatlonal 
stresses, the functional dependence of the matrix material properties on 
temperature and moisture content. In other words, the elastic or plastic 
properties of any matrix material finite element are automatically computed 
to reflect the state of stress and the environmental conditions of temp- 
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orutiiro and humidity. The ndjuatnumt of matorlnl proportion la ir»(}orpor- 
acod with the Incremental loading technique that la employed In thla pro- 
gram. Once the Initial temperature, molature content, and/or elaatlc 
atreaa level for the continuum have boon aped fled, additional loada, be 
they mechanical or environmental, are Introduced In Incromento amall enough 
to permit close approKlmntlon of the nonlinear matrix material propcrtlea 
by amall linear aegmenta. A more detailed description of thla technique 
la presented In Section A. 5 of Appendix A. 

The bulk of the elaatoplaatlc formulation In the present analyala 
program atema from previous work done by Adams [2, 3, 18]. The dovolopmont 
of the hygrothermal loading and material propcrtlea dependence was the 
subject of Miller' a Ph.!). reaearch (1], and the addition of crack Initiation 
and crack propagation capability, to bo dlacuascd In detail In the following 
section, was performed by the present writers. The baala for this program 
la the procedure set down by Zlenklewlci?; |20j, and In fact, the primary 
organisation and flow of the present computer program closely follows the 
suggestions of Appendix A of that text. Thla flow and organisation has 
subsequently been rntlior severely altered to Include crack propagation 
capability . 

The finite element used In this study Is a modified version of the 
familiar constant strain or simplex triangle. For this element, a linear 
dlsplaeoment field within each element Is assumed, to arrive at a function- 
al representation of the potential energy of the system, as referred to 
In Section 3.1 and described In Section A, 2 of Appendix A. The constant 
strain triangular element has some well-known limitations, but for the 
purposes of mlcromechanlcs analyses. It Is an acceptable, economic, and 
powerful tool. The trade-offs Involved In the choice of tlie constant strain 
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triangular oloment Instead of one oC the higher order finite elements la 
covered quite well by Miller in Chapter 3 of Reference [11. A brief 
description of the formulation of the constant strain triangular element 
in its generalized plane strain form is presented in Sections A.l through 
A. 3 of Appendix A. The manner in which loads are introduced to the finite 
element model is described in Section A. 4 of Appendix A. A discussion of 
the concepts and procedure used to analyze the isotropic matrix material 
in the plastic zone, and the temperature/moisture dependence of the matrix 
material properties, is presented in Section A. 5 of Appendix A. In the 
following Section 3.3, the concepts and methodology used to simulate the 
Initiation and propagation of cracks in the composite material are briefly 
outlined, and the analytical relations along with the required sequence 
of operations in the computer program to accomplish this are presented. 


3.3 CRACK INITIATION AND PROPAGATION 

The purpose of the present study was to investigate the effects of 
flaws in unidirectional boron/aluminum composites, with the eventual goal 
of predicting the strength of such composites given a certain statistical 
distribution of Internal flaws. These defects manifest themselves in two 
forms: a discontinuity in one or more boron fibers, or a localized void 

in the aluminum matrix. The loading condition of primary Interest is 
that of tension applied parallel to the fiber axes. With suitable modi- 
fication, a so-called longitudinal model was analyzed with the mlcromechan- 
Ics program in its original form [22]. This permitted modeling of the 
flaw, generally a fiber discontinuity, and an evaluation of the resulting 
localized stress concentration and the local plastic deformation it caused. 
The redistribution of the load to the broken fiber could also be character- 
ized, but only up to the point at which a matrix clement failed (crack 
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initiation). What was neceaaary for further study of the load capability 
of the flawed composite was a crack propagation scheme. This capability 
would also permit a characterization of the energy required to isolate 
the defect in a "zone" of plastically deformed matrix material, or alter- 
natively. the total energy capacity oC the system at the point of cata- 
strophic failure. 

The approach to crack inititation and propagation taken here is known 
as the "failed element" approximation as employed by Adams [2. 3]. l^en 
an element in an area of high stress exhausts its strain energy capacity, 
it fails. From this, we assume that a "crack" has formed and has the 
dimensions of the failed element. This approximation has two implications, 
the most important of which is that a finite amount of material is removed 
from the system, which in an actual material Is not the case. The other 
is that the crack Is not likely to close up on Itself In subsequent load- 
ing because of Its exaggerated width. These effects can be minimized to 
a practical degree by making the finite element grid very fine and uniform 
In the area of anticipated crack initiation. 

It is not enough to simply delete an element from the finite element 
grid when it reaches its ultimate stress. The finite element method in- 
volves the maintenance of force equilibrium at every node point in the 
array, as discussed in Section A. 3 of Appendix A. This equilibrium must 
be maintained when an element falls or unloads. Thus, to represent the 
unloading due to element failure, node point loads which are equal and 
opposite in sense to those equivalent to the state of stress within the 
element at its failure level must be applied at its node points. In 
addition, the failed element's material properties must be set to zero, 
so that the element makes no further contribution to the global stiffness 
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matrix, and all of its computed values of stress and strain are set to 
zero. This insures that the element is completely unloaded and that no 
stresses will be developed in it in subsequent load increments. 

The "reaction loads" applied to the node points of a failing element 
are computed in the following manner. Given the state of stress within 
an element at the time of failure, i.e., cfy, Py, r^y, the statically 
equivalent forces acting at the mid-sides of the element, as illustrated 
In Figure 3.1 can easily be computed: (See for example pp 40-43, Refer- 

ence U91.) 



FIGURE 3.1 

Force-Stress Relationships 

fi’^ = ox^ji + TxyJ^ij 

f/j " <^yHi + "xyyjl 
fjk " °xykj + ^xy’^jk 
fjl = OyXjk + TxyFkj 
file = OxFik ^ T^xyXki 
fik = ^y^ki + "xy^ik 
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where the quantlclcea xij» etc. are aa defined In Figure 3.2. 
which ideneifiea the important geometrical parameters of the generalized 
plane strain triangular element. 



FIGURE 3.2 

Typical Triangular Finite Element of Unit Thickness 

These forces, when translated to the node points with their directions 
reversed, are the reaction loads required for the unloading of the element, 
and are shown below and Illustrated in Figure 3.3. 

■ -W^^j * fi*k) 

k][ - 

Rj . -Wi) * fj-k) 

- -«fi5 + '/k> 

'k ■ + 'jk> 

< ■ -'“'ik ♦ '/k> 
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Node Point Reaction Loads Statically 
Equivalent to Stresses In the Element 

In the present analysis, element failure can occur in one of 
two modes: when both the computed octahedral shear stress and the plas- 
tic octahedral shear strain reach their maximum allowable values (maxi- 
mum dlstortlonal energy criterion) , or when the hydrostatic tensile 
stress In an element exceeds the tensile ultimate strength of the mater- 
ial. This second failure criterion Is also known as failure due to ulti- 
mate cleavage, and failure occurs whenever a tensile principal stress 
exceeds the ultimate tensile strength. 

Although loading Increments are kept small once elements begin to 
enter the plastic region, It is unlikely that an element will fall ex- 
actly at the maximum value of an applied load increment. That is, the 
load Increment will probably be more than sufficient to cause failure 
In the element. For this reason, it Is deslreable to scale down the 
load Increment to the point of element failure. This fraction of the 
load increment required to cause an element to fall Is referred to as 
the "Ratio" and its definition is clarified in Figure 3.4. 
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FIGURE 3.4 

Determination of Increment Scale 
Factor (Ratio) for Element Failure 


Once one or more elements have failed during a given loading incre- 
ment, the crack propagation procedure is set into motion. This procedure 
will be outlined below and is illustrated in the flow chart shown in 
Figure 3.5. 

The first operation to be performed in effecting crack growth is to 
identfy which element failure in a single load increment, if indeed more 
than one element has failed, required the greatest portion of that load in- 
crement to cause failure. The Ratio, as defined in Figure 3.4, is then com- 
puted for this element. This quantity is the fraction of the load incre- 
ment that Is necessary to fall all of the elements in that Increment. It is 
now necessary to multiply all of the element stresses and nodal displacements 
that were computed ns a result of the load Increment application by this max- 
imum value of Ratio. These reduced incremental quantities are then added to 
the sums of the stresses, displacements and strains in the normal manner 
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to give "adjueted" quantities. The nent atop la to compute the nodal 
point reaction loads of the failed elements, using the adjusted accum- 
ulated stress values and the relationships of Equations 3-1 and 3-2. 

These forces become the "adjustment loads" which must be applied to the 
finite element model with the stiffnesses of the failed elements deleted 
from the global stiffness matrix (see Section A. 3 of Appendix A). In 
this manner, the energy required to cause the element failure Is redis- 
tributed Into and absorbed by the remainder of the region being studied. 
This results In new Increments of stress, strain and displacement being 
added to the accumulated values, but the operation Is not counted as a 
load Increment by the program. Of course, this adjustment loading may 
result In the failure of additional elements adjacent to those that have 
already failed, causing the crack to grow without the application of 
additional external loading. This Is similar to the behavior of a crack 
that has reached its critical length. When this occurs, all the elements 
that have failed have nodal reaction loads computed from the stresses re- 
sulting from the adjustment load application, l.e., no ratlolng takes 
place. These elements are then deleted from the global stiffness matrix 
and their nodal reaction loads are applied to the model consisting of the 
remaining elements. This procedure Is repeated until no further element 
failures occur as the result of an adjustment load application, or until 
catastrophic failure of the finite element model occurs. Catastrophic 
failure Is assumed to have occurred when a progression of failed and 
deleted elements results In the division of the model Into two segments. 
This violates the boundary conditions of the analysis scheme and the pro- 
gram terminates, printing the stresses, strains and displacements 
accumulated Just prior to total failure. On the other hand. If successive 
adjustment load application results In no further element failures, the 
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next full load increment is read into the program, added to the accumul- 
ated load, which now reflects n "ratioed" load increment, and the anlysis 
continues as it did before element failure, l.e., if no elements fell, 
another load Increment Is applied, or if failures occur, the entire crack 
propagation procedure outlined above Is repeated. 

3 . 4 DEVELOPMENT OF THK BROKEN FIBER. LONGITUDINAL SFXTrON MODEL 

There are two primary reasons for the development of n longitudinal 
model. One Is to permit study of localized stress concentrations, the 
resulting elastic-plastic behavior of the. aluminum matrix, and subsequent 
crack propagation in the area of boron fiber flaws. Another Is to charac- 
terize the load carrying capability of a flawed boron fiber as a function 
of distance from the location of the fiber flaw. 

These two considerations lead to the most important aspects of de- 
signing the longitudinal model, l.e., geometry, finite clement grid re- 
solution, boundary conditions in the vicinity of a flaw, and spacing of 
the boron fibers in the model. The problem of fiber spticlng will be dis- 
cussed first. 

A typical cross section of a unidirectional, square array, boron/ 
aluminum composite as sho\tfn In the figure below, the section being perpen- 
dicular to the fiber axes. A longitudinal finite element model attempts 
to represent the composite in a plane oriented perpendicular to this sec- 
tion. A longitudinal model of a section parallel to the x or y axes, 
through the centers of the fibers, would be representative of the minimum 
distance between fibers. A section cut at 45“ to the x-axis and through 
the fiber centers would depict a maximum fiber spacing situation, lilhen 
one of these fibers is broken, the load it carries decays to zero at the 
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broken surface, assuming that the boron-aluminum interface remains intact. 
At the flaw site, the fibers adjacent to the broken fiber, and to some ex- 
tent the surrounding aluminum matrix, must absorb the load that the broken 
fiber would have otherwise carried. The aluminum transfers this excess 
load back into the broken fiber via a shear mechanism so that at some dis- 
tance from the fiber break, that fiber is again fully effective in carry- 
ing load. It is logical to presume that the amount of aluminum between 
the boron fibers will have an effect on this load transfer mechanism. To 
characterize the effects of variation in fiber spacing, two longitudinal 
models were studied, one representing a 90*’ section cut of the transverse 
cross section, and another representing a 45° section cut. These two 
models are shown in Figures 3.6 and 3.7. Note that the fiber diameter 
dimensions have been normalized to unity. In Figure 3.7 the effect of 
the 90° section cut in diminishing the amount of local aluminum matrix is 
shown quite clearly. The size and aspect ratios of the fiber elements are 
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exactly the same as those of Figure 3.6, but the aluminum elements of the 
90** section model are so compressed that the element numbers, which are 
identical to those of Figure 3.6 have been eliminated for clarity. 


The second problem that must be solved in the finite element modeling 
of a broken fiber in a composite is the geometry in the area of the fiber 
discontinuity-matrix interface. The efforts of the present investigators 
to resolve this problem have been evolutionary in nature and were described 
in detail in Section 5.2 of Reference [23]. 

The results of this evolutionary process are the finite element models 
shown in Figures 3.6 and 3.7. Some of the features of this model can be 
Illustrated by referring to Figure 3.8, a section of the model represent- 
ing the discontinuous fiber and the local aluminum matrix. Note that the 
free surface of the discontinuous boron fiber is extended out to the model's 
plane of symmetry by the addition of fiber element number 157. Longitud- 
inal models run without this element developed extremely high stress 
levels in aluminum element number 1 at very low levels of applied load. 

While this situation might occur in a boron /aluminum composite if the 
fiber were discontinuous and its ends separated by some finite amount at 
the time of fabrication, it was felt that the present configuration, as 
shown in Figure 3.8, was of more general use. This configuration might 
be used to represent a fiber that has broken during fabrication, or, if 
flawed locally, failed at a very low loading of the unidirectional com- 
posite. Note the small size and uniformity of the aluminum elements in 
the area of the end of the discontinuous fiber. This is to permit a closer 
approximation of actual crack growth, as discussed in Section 3.3. In 
order to prevent the broken fiber from having any load capability at x ° 0, 
node point 12 has been relased from its x-dlrectlon fixity. This results 
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In a free surface for element 2 along the y-axis. or in other words, an 
Initial crack in the aluminum as might actually occur from the initial 
release of energy caused by the fiber failing at a low load level. 

3.5 DEVELOPMENT OF THE TRANSVERSE SECTION MODEL 

Finite element modeling of a transverse section of a unidirectional 
boron/aluminum composite is fairly straightforward. However, the need to 
study the influence of a reduced load capacity in one fiber on its neigh' 
boring fibers requires that a minimum section model such as that shown in 
Figure 3.9 be employed. This model represents the first quadrant of a re- 
peating unit cell of a rectangular array of fibers. In effect, a flawed 
fiber can be considered to be surrounded by eight other fibers in the 
array. A model of this type can easily lead to a great number of elements, 
and attempting to Increase the resolution of the grid at selected locations 
often results in a very large bandwidth of the overall stiffness matrix 
for the finite element model. The transverse model developed and reported 
in the First Quarterly Report [22], proved to be too large for the University 
of Wyoming's present computer. This resulted in the development of the 
model shown in Figure 3.9. Also, during the last year, the mlcromechanlcs 
analysis computer program has been modified to permit the inclusion of as 
many as four different materials, each having different properties. This 
capability is essential in representing breaks or flaws in the boron fibers 
of a transverse model via reduced stiffness. 
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SECTION 4 


MATERIAL PROPERTIES 

In modeling the boron/aluminum composite, the boron fibers have been 
treated as brittle, linearly elastic materials with isotropic strength 
and stiffness properties. The aluminum matrix has also been considered 
to be isotropic, but Is modeled as an elastoplastic material. To accomp- 
lish this, the actual stress-strain curve of the aluminum alloy selected 
is input to the analysis by curve fitting via a Richard-Blacklock two- 
parameter equation, as discussed in Section A. 5. Thus, at any load level 
the tangent modulus for any given element can be computed. This makes 
possible an accurate representation of the plastic deformation of the 
matrix. 

Although the nonlinear material properties of any matrix material, 
e.g., another aluminum alloy, can readily be incorporated in the analyses, 
a 6061-T6 aluminum alloy at 75°F was chosen for the Initial studies. The 
material properties shown in Table A.l were obtained from Reference [24]; 
the full range stress-strain curve for determining the curve fit para- 
meters used is shown in Figure 4.1. 


TABLE 4.1 

Aluminum Matrix Material Properties - 6061-T6 Alloy [24] 


E ® 10.0 X 106 psi 
V “ 0.33 
Fty = 36000 psi 
F^u - 45000 psi 
a “ 13.0 X 10-6 in./in./“F 


Young's Modulus 
Poisson's Ratio 
Tensile Yield Strength 
Tensile Ultimate Strength 
Coefficient of Thermal Expansion 
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Stress, 1000 psi 


Th 0 boron fiber properties Indicated in Table 4.2 were obtained from 
Reference (25]. 


TABLE 4.2 


Boron Fiber Moterlol Properties [25] 

Young's Modulus E ■ 60.5 x 10^ pal 

Poisson's Ratio v ■ 0.130 

Tensile Strength ■ F^X - 500,000 pel 

Coefficient of Thermal Expansion a ■ 9.0 x 10"^ In./ln./^F 
Ultimate Strain - 8264 x 10"6 In. /In. 



Strain, in. /in. 

FIGURE 4.1 

Typical Full Range Stress-Strain Curve For 
6061-T6 Aluminum Alloy at Room Temperature [24] 
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SECTION S 


NUMERICAL RESULTS 


Tho Introduction of crack Initiation and growth capability to tha 
micromechanics program has made possible a much more extensive study cf 
the effects of o discontinuous fiber in a boron /aluminum composite. How- 
aver, a numerical difficulty was revealed in the crack propagation pro- 
cedure shortly before the end of this report period. Debugging of the 
subsequent corrected computer progr^ has prevented the presentation of 
the crock growth results in this report; they will be presented ii\ full 
detail in the next quarterly progress report. In spite of this difficulty, 
the results of loading both the AS” section and the 90” section models to 
the point of first element failure are available, and very informative 
significant differences in the load carrying capability of the 43” section 
model versus the 90” section model hove been revealed. It is antici- 
pated that future analyses with crack propagation will generate even more 
important insights concerning the generalized plane strain treatment of 
this problem. 

The 45” section model and the 90” section m.odel, as described in 
Section 3.4, represent a unidirectional composite with 25 percent of the 
fibers containing a break. Larger models with more unbroken fibers can 
easily be studied If desired. Each of these models has been loaded to 
the point of crack initiation (first element failure), and the difference 
In load capability of each at this point is considerable. In both cases. 
Initial plastic deformation occurs at very low load levels, and due to 
considerable plastic strain capability of the aluminum matrix, as shown 
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In Plgura 4.1, eho loading Incromonca bayond choc point era kape vary 
asall, i.at, 1000 to 3000 pai. Conaidarabla locolizad plaatle daforaation 
waa devalopad in tha ragion of tha flbar discontinuity befora an alamant 
failura occurrad. Tiiis in turn ganarally triggarad tha failura of aavarol 
mora inatriK alamanta bafora anothar loading incramont could bo applied. 

In the paragraphe that follow* the roaulta of loading the 45" aection and 
tha 90" aection longitudinal modela up to the point of crack initiation 
will be studied. In Section 5.3* the results of axial loading of the 
transverse model and simulation of a broken fiber are reviewed. 

5.1 AXIAL LOADING OF THE 45" SECTION LONGITUDINAL MODEL 

The 45" section longitudinal model* as illustroted in Figure 3.6* 
was loaded in the direction of the fiber axes (x-direction) to a level 
of 66*076 pel average applied stress* at which point element number 2 
failed. Initial plastic deformation was also observed in element number 
2, at an average applied stress level of 17*000 psl. Plastic deformation, 
in this analysis, is defined as the point at which the octahedral shear 
stress-octahedral shear strain curve becomes nonlinear. For the 6061-T6 
aluminum alloy used in this study, the octahedral shear yield strength 
is 16*970 psl. Element number 2 would be the first to experience plastic 
deformation, for the reasons cited in Section 3.3. 

Loading was increased monotonlcally until element number 2 failed* 
at an average applied stress level of 66,076 psl. At this level there 
was considerable plastic deformation in the region of the fiber discontin- 
uity, with several elements adjacent to number 2 very near failure. By 
using the plotting capability of the micromechanics computer program, this 
zone of plastic deformation at the point of imminent failure of element 
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nuiBbor 2 can be soon. Xn Figure 5.1, contours of conacant occahodral 
■haor atreas for che aluminum matrlK have bean plotted on an outline of 
the 45* aection finite elemn'**' model. Note that the contour valuea are 
for octahedral ehear atreaaea normalized with reapect to the yield atrength, 
16,970 pal. Thia meana that any areo encloeed by contour linea that are 
greater than or equal to one haa been plaetically deformed. Aa Figure 5.1 
ahowe, conalderable plaatic deformation haa occurred at the 66,000 pai 
applied atreaa level. In addition, almilar contour plota for octohedral 
ahear atrain, maximum principal atreaa, minimum principal atreaa, and in- 
plane ahear atreaa (Tjty) ore provided in Flgurea 5.2, 5.3, S.4» and 5.5, 
respectively. These Illustrate very well the general state of stresa in 
the aluminum matrix at the load level of imminent crack initiation. 

A study was also made of the loading intensities in the boron fibers 
at various load levels, in an attempt to further understand the mechanisms 
of load transfer to the discontinuous fiber as plastic deformation increases. 
In Figures 5.6 and 5.7, the relative loading intensities of the three 
boron fibers in the model are compared as a function of the axial distance 
from the flaw site. The two lood levels chosen are 16,000 psi for Figure 
5.6 and 66,000 psi for Figure 5.7. From these figures it can be seen 
that the end of the broken boron fiber carries no axial stress, as it 
must, and that the load level in the adjacent fiber rises abruptly in the 
vicinity of the fiber break. Of particular Interest is the relative 
loading of the discontinuous fiber as the load level increases. It will 
be noted that the slope of the load curve for the discontinuous fiber at 
che flaw site decreases as che applied load Increases. In addition, at an 
applied stress level of 16,000 psi, Che broken fiber has attained a stress 
level 78 percent that of a remote fiber at 4.55 fiber diameters from the 
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Octahedral Shear Strain Contours Just Prior to Crack Initiation, 45“ Section 

Longitudinal Model, ksi * 





Contours of Mininiuin Principal St 







DISTANCE IN FIBER DIAMETERS FROM THE PLANE OF FIBER BREAK 

FIGURE 5.7 

Fiber Loading Versus Distance From the Broken Fiber Surface, 45® Section Hodel, 

66,000 psi Average Applied Stress 




flaw site, while at an average applied stress level of 66.000 psl. the 
broken fiber has been loaded to only 73 percent at 4.55 fiber dian^ters. 
These values imply that at a aufficlantly large distance from a flaw 
site, all the boron fibers will be equally loaded, but that the distance 
required for this condition to exist increases as the average stress level 
increases. Crack propagation will no doubt cause the distance from the 
crack at which the broken fiber is again fully loaded to be much greater 
still. 

As an illustration of the energy absorption capacity of the 45“ 
section broken fiber longitudinal model, the composite stress versus 
composite strain has been plotted up to the point of crack initiation, as 
shown in Figure 5.8. It will be noted that the curve is essentially 
linear, in spite of the considerable plastic deformation shown in Figure 
4.1. This is due to the fact that the boron fibers, whose axial modulus 
is very much higher than that of the aluminum, are essentially linear in 
their stress-strain response. In fact, even at the plane of fiber dis- 
continuity, the boron fibers are carrying 76 percent of the applied load. 

At two fiber diameters fro'o that plane, 90 percent of the applied load is 
being carried by the boron fibers. The overall stiffness of the 45® 
section model, i.e., the shape of the stress-strain curve as sho\im in 
Figure 5.8, proved to be 36.0 x 10^ psl, or about 60 percent of the axial 
stiffness of the boron fibers. 

5.2 AXIAL LOADING OF THE 90® SECTION LONGITUDINAL MODEL 

In loading the 90® section longitudinal model, Figure 3.7, parallel 
to the boron fiber axes, plastic deformation was first observed in element 
number 2 at an applied average stress level of just under 10,000 psi. This 
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Compotitt ttroin, in./in.) 


FIGURE 5.8 

Plot of Applied Stress Versus Composite Strain 
for the 45** Section Longitudinal Model 

is substantially lower than the 16,000 psi level that was necessary to 
cause plastic behavior in the 45** section model, and is due primarily to 
the fact that the apparent volume of aluminum matrix available to trans- 
mit load between the discontinuous and Intact boron fibers is only about 
30 percent of that seen by the 45** section model. 

Loading of the 90** section model was Increased until element number 
2 exhausted its strain energy capacity at 40,150 psl applied stress. Again, 
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HtresB and strain contour plots were penorated for the load level which 
was reached immediately prior to the first failure, and these are shown 
in Figures 5.9, 5.10, 5.11, 5.12, and 5.13 for the normalized octahedral 
shear stress, octahedral shear strain, maximum and minimum principal 
stresses, and the in-plane shear stress, respectively. It will be noted 
that the zone of plastic deformation, as defined by the 1.03 contour in 
Figure 5.9, is considerably less extensive than was observed in the 45° 
section model. Figure 5.1. From preliminary results of crack propagation 
computer runs, it is quite apparent that the pattern of deformation, flaw 
growth, and fiber loading differs quite significantly between the 45° and 
the 90° Section longitudinal models. Further, the differences in load 
levels necessary to initiate plastic deformation and later, crack formation, 
for the two models Indicates that the proximity of adjacent fibers to a 
flawed or discontinuous fiber has an Important effect on the behavior of 
the total composite material. This problem will be discussed further in 
Section 6. 

Another difference in the response of the 90“ section model when 
compared with the 45° section model Is shown in Figure 5.14, in which the 
composite stress-strain response of the 90° section model is plotted up 
to the point of crack initiation. The composite axial modulus of the 90° 
section model is found to be 49.0 x 10^ psi, using Figure 5.14. This is 
significantly stiffer than the 45° section model, as would be expected 
in view of the larger boron fiber volume fraction of the 90° section model. 
The strain energy necessary to initiate failure in the 90° section model 
was determined to be 36.1 in. -lb. /in., while in the 45° section model, 

123.2 in. -lb. /in. of energy was absorbed before crack initiation. 
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Longitudinal Model (10“3 in. /in. 





Longitudinal Model (ksi) 
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FIGTE 5.14 

Plot of Applied Stress Versus Composite 
Strain for the 90° Section Longitudinal Model 


5.3 AXIAL LOADING OF THE TRANSVERSE SECTION MODEL 


The transverse section model, as Illustrated in Figure 3.9, was loaded 
In the out-of-plane or z-dlrectlon, both In the unflawed condition and with 
one broken fiber. In order to approximate the broken fiber plane of the 
longitudinal models, the transverse model was run with the stiffness prop- 
erties of one of the fibers reduced to zero, as shown In the following sketch 



The percentage of the total load carried by the fibers of the un- 
flawed transverse model proved to be 89.0 percent t which compares very 
well with the percentage of load carried by the fibers of the 45“ section 
model, which was 89.1 percent. When the transverse model was run with 
one fiber deleted, it was found that the three remaining fibers carried 
85 percent of the applied load, which is in good agreement with the 90“ 
section longitudinal model, which showed 86 percent of the load in the 
boron fibers at the plane of discontinuity. 

The stress distribution in the aluminum matrix for the transverse 
section model is shown in Figures 5.15 and 5.16. These plots represent 
the stress in the fiber axis direction (Oj. in this case) for the matrix 
material along the y-axis and along a line 45“ to the y-axis. Both the 
stresses for an unflawed composite arid one in which the fiber centered 
at X « y - 0 has an effective modulus of zero are shox^m. The applied 
stress is 2,000 psl. As these plots show, there Is a definite change 
in stress levels when one fiber is deleted; the stress gradients between 
the fibers are quite low when compared to those seen in the longitudinal 
models. This of course is due to the fact that the transverse model can 
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Transverse model 
with one broken fiber 


Broken fiber 
boundory 

Unflowed tronsverse 
model 


0 0.05 O.iO 0.15 0.195 

Distance between Fibers, (fiber diameters) 


nOlIKli 5.15 

Axial Stress Distributions in the Aluminum Matrix, Transverse 
Model, 90" Section (minimum distance between fibers). 

Applied Stress a.. = 29,000 psi 


Axial Stress in the Aluminum Matrix, G^(IO^psi) 


Broken fiber 
boundory 


Transverse model with 
one broken fiber 


Unflowed transverse 


0.0 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 

Distance between Fibers, (fiber diameters) 

FIGURE 5.16 

Axial Stress Distributions In the Aluminum Matrix, Transverse 
Model, 45” Section (maximum distance hetween fibers), Applied 
Stress Oj 5 = 29,000 psl 
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only reproBont dll;£ercncc8 In material properties caused by broken or 
reduced capacity fiber elements, but cannot account for the physical dla- 
contlnuity of a fiber and the resulting stress concentrations. In fact, 
for both the unflawcd case and the model In which one boron fiber Is de- 
leted, the variation In stresses In the axial or z-dlrectlon Is no more 
than about 10 percent throughout the aluminum matrix, and considerably 
less then this in the boron elements. 

The In-plane stresses, l.e., o^, Oy, and Txyt show considerable vari- 
ation In the aluminum matrix, but arc about 10 orders of magnitude less 
than the axial stress. Contour plots for constant values of In-plane 
shear stress, and maximum and minimum principal stresses, arc shown In 
Figures 5.17, 5.18, 5.19, 5.20, 5.21 and 5.22, for both the unflawed and 
flawed longitudinal models, to show the Influence of deleting one fiber. 
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FIGURE 5.19 

Shear Stress Contours (ksl), Unflawed Transverse Section Model, Applied 

Stress Ojj ■ 32 ksl 








FIGURE 5.22 


Shear Stress Contours (ksi), Transverse Section Model With One Broken 

Fiber, Applied Stress Sg • 32 ksi 





SECTION 6 


FUTURE WORK 


Debugging of the updated and corrected crack propagation version of 
the mlcromechanlcs computer program is nearly complete. Production runa 
using this program will soon be performed on the 45° and the 90° section 
models. Behavior in terms of crack propagation and energy absorption will 
be determined, and attempts to correlate these results to the experimental 
results of other investigators will be made. From experience to date with 
crack propagation analysis schemes, a further refinement of the 90° section 
longitudinal model appears desirable. Specifically, the aspect ratios of 
some of the finite elements are presently rather large; the region of uni- 
form and small elements required to more accurately study crack growth will 
be extended. Once the energy absorption characteristics of these longi- 
tudinal models have been determined, several additional areas of interest 
will be explored. 

Transverse loading of both longitudinal and transverse section models 
containing flaws will be of Interest, particularly If these loading con- 
ditions are preceded by temperature gradients and hydrostatic loading 
Increments to simulate the fabrication process for boron/aluminum. This 
"preconditioning" provides a prediction of the residual stresses and 
strains In the composite prior to mechanical loading, and could have a 
significant effect on the overall Influence of Internal flaws. 

Finally, it has not been possible, to this point, to relate the 
stress concentration effects seen in a longlti^dlnal model to a transverse 
section model. This would not be a serious problem were It not for the 
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illffercmcsB In behavior obaervod In axial loading of Cho 90* and Che 45* 
aectlon modela. Purthertnore , in examining a tranaverae aectlon of a typical 
square array unidirectional compoaite» aa llluatrated in Section 3>3» one 
can aee that the flberi^ in cloaeat proximity to a broken fiber, e,g., 
those represented by a 90* section model, will be affected in turn by 
their closest unflawed neighbors, some of which ore the continuous fibers 
represented in a 45* section mod«>l. In effect, all of the fibers surround- 
ing a broken fiber are coupled in the mechanism of load redistribution, 
and crack propagation is likely to be too complex to be characterized by a 
two-dimensional formulation. Until more extensive experimental data ore 
available to verify and improve the present two-dimensional formulation, 
a three-dimensional finite element micromechanics analysis appears to be 
the logical approach to further understanding of energy absorption in 
metal-matrix composite'/’.. Such a program is in the final stages of devel- 
opment at the University of Wyoming, and during the second-year study, pre- 
liminary investigations of its applicability to the present problem will 
be made. 
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APPENDIX A 


FINITE ELEMENT FORMULATION OF THE 
MICROMECHANICS ANALYSIS COMPUTER PROGRAM 


A.l COORDINATES AND BOUNDARY CONDITIONS 

To aid in the discussion of topics to follow, Figure A.l is presented 
as representative of a typical repeating unit of a composite material 
under analysis, and Figure 3.2 defines some of the geometrical parameters 
describing a typical triangular element. 



AT z = 1 , 8z = constant 
Tz= Tyz =0 

FIGURE A.l 

Configuration of the Area of Interest for a Typical 
Mlcromechanlcs Analysis, Including Boundary Conditions 


As discusHCii by Miller (pp. 32-42 of Reference flj), the fiber 
packing geometry Is assumed to be symmetric, and the repeating first 
quadrant unit of the composite continuum represented by Figure A-1 must 
remain rectangular when any combination of normal tractions, Ox* 
and Oz are applied to the three orthogonal planes. This results In the 
specification of unique boundary conditions to maintain the symmetry and 
continuity of the material under Investigation. These boundary conditions 
arc explained In detail by Miller and Adams [1], and by Adams [2J^ and 
sumnuirlzed In Figure A.l. 

A. 2 GENERALIZED PLANK STRAIN 

In Miller's fl] formulation of the governing constitutive equations 
for the finite elements used In this analysis, generalised plane strain 
conditions were assumed In order to reduce the analysis to a quasi-three- 
dlmenslonal problem. In past micromechanics Investigations, finite ele- 
ment models representing a transverse section of a unidirectional com- 
posite were treated as ordinary plane strain problems. Under these con- 
ditions, the body under consideration is assumed Lo be thick, l.e., the 
axial dimension Is much larger than the transverse dimensions, and the 
displacement In tlint direction is assumed to be zero, resulting In zero 
strain In that direction, l.e., 

Ez “ Yxz= Yyz “ 0 (A-l) 

For isotropic materials, this condition permits an Induced normal stress 
In the axial direction of 

Oy. = V(0^ + Oy) (A-2) 

Miller's primary reason for incorporating generalized plane strain con- 
ditions was to allow axial loading, l.e., loading parallel to the fiber 
axes, the z-directlon in his transverse section finite element grids (see 
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Figure A-1) . In addition, when a composite is loaded in the transverse 
direction with no longitudinal tractions being applied, the generalized 
plane strain conditions allow a normal strain to develop along the longi- 
tudinal axis, resulting in zero stress in the longitudinal direction. 

Lekhnltskli [26j defines generalized plane strain in a very general 
manner, stating that all of the strains associated with the z-axls di- 
rection can be nonzero constants. Including the shear strains, Yxz 
V y 2 « In the current analysis, a somewhat less general form of generalized 
plane strain has been employed in which it is assumed that only the normal 
strain, is nonzero, i.e., Yxz nnd Yy^ n^e assumed to be zero. This 
definition requires that all planes perpendicular to the z-axls direction 
be a linear function of the position coordinates in that direction, i.e., 

(i> “ Kz (A-3) 

where k is a constant. In other words, some constant strain, exists 
in the z-axls direction 


3(0 

3z 


= K 


(A-A) 


This means that the displacements of all of the node points in the z-axis 
direction are identical. It is also Important to note that the axial 
stress, is uncoupled from the transverse stresses, Oy, and Tj^y. 
Keeping this in mind and referring to Figure 3.2, the following points 
are noteworthy: 

a) Node points 1', j', and k' are required to have the same dis- 
placements in the x and y directions as points i, j, and k. 

b) Node points 1, j , and k have Identical displacements in the 
z-dlrectlons, while from symmetry considerations, i' , j', and k’ 
have zero displacements in the z-directlon. 

c) Since each element is assumed to be of unit thickness, the dis- 
placements of the unprimed node points in that direction represent 
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the strain, 

A. 3 FORMULATION OF ELEMENTAL STIFFNESS MATRICES BASED ON GENERALIZED 

PLANE STRAIN CONDITIONS 

The full development of the basic finite element method is avail- 
able in complete detail in several texts (References [19-21], for example), 
and a complete examination of the formulation and assembly process re- 
sulting in a set of simultaneous algebraic equations will therefore not 
be included here. Only a brief overview of the process is presented so 
that the unique capabll Itles of the present analysis scheme might be better 
explained. 

Once the assumed displacement field within each element has been e.:- 
pressed in terms of unknown node point displacements (see Figure 2.1), the 
strain at any point within the element can be expressed as 

Hi * (A-5) 

with i = 1,2 

Uj^ = u = X- displacement 
U 2 ■ V = y- displacement 

XI = X 

X 2 = y 

As a result, a general expression relating the strains to the displacements 
in any element, i, can be written as 

{G)i = [Bli (A-6) 

The matrix is a set of geometric parameters relating the vector of 

the node point displacements, {iS)i, to the strains in the i-th element. 

The form of [8]^, generally knov\m as the "shape" matrix, is dependent on 
the form of the assumed displacement field. 

The stresses and strains are related by an appropriate constitutive 
relationship, Hooke's law being a familiar example. The general Duhamel- 
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Neumann form of these constitutive relationships (20, 27, 28] is 

(o) - [H] ((c) - {cq}) (A-7) 

Wlicre (o) and (c) are the stress and strain vectors, respectively, and 
[H] is a matrix containing tht\ appropriate moterial properties relating 
the two. The elements of the initial strain vector, {c^}, are dilata- 
tional strains induced by thermal or moisture changes. This topic will 
be discussed in greater detail In Section A. 5. 2 of this Appendix, but 
for now it is noted that the temperature and moisture sensitive matrix mat- 
erials have an order of symmetry which Is at least orthotropic, so that 
dllatatlonal shear strains cannot exist. That is, the initial strain 
vector can be expressed as 


ieo> ■ 


r.. "I 





(A-8) 


As described by Zlenkiewicz [20] for instance, an expression for the 
strain energy within an element can be written using the relationships 
given above, leading to an expression for a potential energy functional 
In terms of strain energy and forces acting at the node points. Mini- 
mization of this functional results In a set of simultaneous algebraic 
equations relating node point forces to node point displacements. 
Typically, 


- [k]^[S}^ + (F^ )i (A-9) 

o 

where {F)i is the vector of forces acting on the node points of element 

1, and [k] is the stiffness matrix for element 1. The elements of the 

vector {F are the forces acting on the node points resulting from 
^0 

initial dllatatlonal strains developed within the element. For constant 
strain elements, the form of [kjj^is 
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(A-10) 


Ikli - Ull (HJi iBJi tiAi 
where Is the thickness of the element and is the area of element i 
in the x-y plane. 

All of the forces on each node point in an array due to each element 
sharing that node point must be in equilibrium. In addition, the summation 
of all such node point forces must be in equilibrium with all externally 
applied loads and specified boundary restraints. This equilibrium require- 
ment and the summation process it entails leads to the formation of a 
"global" stiffness matrix, Ik1» for the entire region under analysis. 

Wliat this essentiall'y does is to assemble all of the equations of the form 
of (A-9) for all elements, resulting in a total set of simultaneous equa- 
tions for the entire area being analyzed. 

{F)« [K] { 6 } + {F } (A-11) 

^o 

For our constant strain element, the displacement fields under gener- 
alized plane strain conditions are 

u * a^ + 32 X + 33 y 

V “ bj^ + b2X + b3y (A-12) 

w ® tez 

where k is a constant. 

For the triangular element shown in Figure 3.2, the coefficients Aj, A 2 , 

A 3 , b^, b 2 , and bj can all be expressed in terms of the node point dis- 
placements in the x-y plane. This leads to the shape matrix as de- 

veloped by Heubner [19], but with a fourth row and a seventh column added 
to express the condition of generalized plane strain: 
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where • a'AT + B'AM 

For a complete dlscusaion of the derivation of the material properties 
matrices shown above, the interested reader will find a thorough and 
easily understood development by Miller and Adams (pp. 47-53 and Appendices 
A and C of Reference [1]). It is to be noted that while these expressions 
for involve the region of elastic or linear stress-strain behavior, 
the values of the material constants E, v, a, 0, etc., can be functions 
of temperature and moisture. 


A4.0 LOAD APPLICATION 


Loads introduced to the finite element array in the micromechanics 
analysis can be In the form of applied mechanical tractions, or arise from 
thermally or moisture induced dilatations. Mechanical loading can consist 
of average applied normal stresses in the x, y, and z directions, as de- 
fined by Figure A-1, for each load Increment, while thermal and moisture 
gradients can be applied at any Increment to reflect environmental changes 
the composite may be subjected to. Mechanical loading will be discussed 
first. 

The application of mechanical tractions to the finite element model 
is considerably simplified by taking advantage of the boundary conditions, 
as specified in Figure A-1, which permits a rearrangement of the global 
stiffness matrix [K], and the total force vector-, {F}, by a method intro- 
duced by Branca [30]. The displacement boundary conditions for the re- ' 
peatlng unit finite element model were specified in order to maintain 
continuity of the material continuum under investigation. Specifically, 
referring again to Figure 2.1, displacements in the x-direction of node 
points along the right-hand vertical boundary must be uniform. Dlsplace- 
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meats in the y-dlrectlon of the upper horizontal boundary roust be uniform, 
and the displacements of all node points in the z-axis direction roust be 
uniform. When the overall force-displacement equation of the system is 
considered, i.e., 

{F} - fKl {«} (A-18) 

one can see that all of the boundary node points Involved in mechanical 
loading will have Identical displacements with respect to the direction 
of the load application. These identical displacements allow combining of 
certain terms in the global stiffness matrix that result In the replacement 
of the applied forces on boundary nodes by zeroes, in the manner described 
by Branca [30]. Successive modification of the global stiffness mat- 
rix for each boundary node point displacement results in the following 
form of Equation (A-18) for the simultaneous application of uniform values 
of Ox, Oy, and for an array of n nodal points 

•^ll *<12 •••kl(2n+23] 

k22 k23 


^0'' 
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> - 








^33 



^«1 ^ 


«2 
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« 2 n 


^ 2 n+l 


»*52n+2»/ 


(A-19) 


where Fx, Fy, and F 2 are the total applied loads In the x, y, and z direc- 
tions, and are defined, for a unit thickness model, as 

Fx “ ^xb 


Fv “ o„a 


(A-20) 


Fz = Ogab 
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where a and b are as defined In Figure 2.1. On pp. 63*67 of Reference [1], 
Miller clearly illustrates, by means of a simple two-element model, the 
manner In which the elements of the global stiffness matrix are manipulated 
to obtain Equation (A-19). The Important point Is that this procedure 
allows a rather straightforward simultaneous application of external trac- 
tions In the directions of the three coordinate axes. 

Thermal-dilatatlonal and molsture-dllatatlonal effects are also In- 
cluded within a given loading Increment. These effects appear as the vec- 
tor {Fp^)^ In Equation (A-9). The elements of this vector represent the 
forces on the node points of an element which arc the result of dllatational 
strains, e^, due to thermal and moisture gradients, as given by Equation 
(A-15). The magnitudes of these Induced nodal forces are a function of 
the shape of the element, [B], and Its materlol properties, [H], l.e,, 

(F^ >1 - (A-21) 

o 

tiTlien a loading Increment Includes a change In temperature or moisture 
content, the node point forces given by Equation (A-21) are calculated 
and then moved to the right-hand side of Equation (A-9) 

(F-F^ (A-22) 

This form la retained when the elemental equations are assembled in- 
to the total global form, as given by Equation (A-11) , so that a single 
loading vector exists for each increment and only one inversion of the 
total stiffness matrix, [K], Is required. 

A. 5 NONLINEAR MATERIAL RESPONSE 

The University of Wyoming micromechanics analysis program models the 
response of the Isotropic matrix of composite materials to both external 
tractions and changes In temperature and/or moisture concentrations. Be- 
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cause the vaet najority of matriv materiale In uac today are capable of 
considerable plostlc deformation before failuroi an elastoplastic analysis 
scheme for deformation due to mechanical loading is ossential for micro* 
mechanics studies. In addition, the effects of temperature and moisture 
(hygrothermal) changes must be considered. Internal stresses Induced by 
the high temperature fabrication and bonding processes that most composites 
are subject to can actually cause material failures before any mechanical 
load is applied, while tne effects of moisture absorption by polymeric 
resins has become an area of major concern for analysts and designers in 
recent years. Past micromechaiilcs studies (2, 3, 7-18] have assumed that 
the material properties remain constant for all states of temperature and 
moisture to which the constituent materials are exposed. This is not the 
case, as these changes in temperature and moisture content not only induce 
stresses, but at the same time significantly alter the material properties. 
This Is also a type of nonlinear behavior and it must be accounted for. A 
brief description of the manner in which the effects of plasticity and the 
hygrothermal dependence of the material properties are Incorporated Into 
the finite element stiffness matrix follows In the next two subsections. 


A.5.1 ISOTROPIC MATERIALS IN THE PLASTIC REGION UNDER GENERALIZED PUNE 
STRAIN CONDITIONS 

Of the many elastoplastic finite element analyses of undlrectlonal 
composites performed to date [2, 3, 8-11, 13, 16, 18], the present proced- 
ure is most like that described by Adams [13]. It 1st modified to account 
for the generalized plane strain assumptions and for the inclusion of the 
effects of environmental changes. Like most methods of approximating non- 
linear material behavior, the present procedure requires that loading be 
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applied In BBall Incremente. Thia nececslcatea the accumulation of dla- 
placecMnta, acreaaea and attains from each succeeding Increment. Xt will 
be shown that this Is possible due to the fact that the material properties 
used In Incremental techniques are linear and therefore superposition prin- 
ciples apply. 

Two basic techniques exist for accounting for the plastic portion of 
each loading Increment; the tangent modulus method and the method of In- 
itial strains. The method of initltal strains uses the elastic material 
properties, as employed In Equation (A-14), throughout the entire loading 
sequence. It accounts for plastic strain by adding an Initial strain to 
the strain vector of each element and then Iterating until equilibrium con- 
ditions are satisfied. The advantage of this method is that the global 
stiffness matrix [K] need only be formed and inverted once. The tangent 
modulus method uses the tangent modulus of the material at a particular 
stress-strain state to define the stiffness of eacn finite element for the 
next loading Increment. While no Iteration is required, a new global 
stiffness matrix must be assembled and Inverted for every loading incre- 
ment. A more detailed description of both of these methods is given by 
Adams (pp. 36-39 of Reference [2]), wherein the tangent modulus method 
emerges as being preferred for materials that exhibit relatively "flat" 
stress-plastic strain curves, which Is the case for many metal matrix ma- 
terials used In advanced composite materials. The current analysis uses 
the tangent modulus method for this reason. In addition, since the mater- 
ial properties must be modified at the beginning of each loading increment 
to account for environmental factors, which also requires the assembly of 
a new global stiffness matrix, use of the tangent modulus method Imposes 
no penalty. 


For each loading Increment beyond the elastic limit, the correspon- 
ing Increment of strain fo' any element can be separated Into an elastic 


(recoverable) and a plastic (irrecoverable) part, or in Adams' notation [2], 

(< 

IJ 

(o) 


eij ' (A-23) 


where the elastic portion, behaves according to the generalized 

Hooke's law. 


(e). l-2v 


1+v 


(A-24) 


Hj'"'- ^kk^lj + ^ ^ij 
where is the Kronecker delta, and s^j is the devlatorlc component of 

the rate of stress tensor, 

1 


Hi - *lj - 3 ^kk^ij (A-25) 

For this analysis, the plastic portion of the deformation is assumed to 
follow the Prandtl-Reuss flow rule [31], 

is,. (A-26) 


Hj - ''“ij 

where X is a positive scalar function. In other words, at any Instant, 
the rate of change of the plastic strain is proportional to the devia- 
toric stress only. The mean normal, or dllatatlonal, strain makes no 
contribution to plastic deformation. Adams (pp. 13-15 of Reference [2]), 
goes through a detailed explanation of the procedure involved in obtaining 
a convenient form of : X 


X = 


2toMt 


(A-27) 


where Tq is the octahedral shear stress 


Tq = (1/3 sijsij)*5 (A-28) 

To is the octahedral shear stress rate of change, and is the tangent 
modulus of the octahedral shear stress-octahedral plastic shear strain 
curve 

2M, = (A-29) 
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a linear relationship. The above equations lead to an incremental con- 
stitutive equation for ><^la8toplaatic material behavior 

(A-30) 


. l+v . V . . *’ij®kl‘^kl 

"ij " If ^Ij - E^’kk^iJ + 

OToMp 


Miller [1) adds an additional term to this expression to account for the 
dilatationul strain increments induced by the environmental changes. 
Thus, 

(A-31) 


• l+v • V • j j. y. i * a: 

‘'ll “ ¥ ‘^ll ■ E '^kk'^ij + + to^ij 

6tq Mt 

Equation (A-31) must be inverted in order to obtain the constitutive re- 
lationship in the form of Equation (A-6) . Details of this inversion pro- 
cedure under generalized plane strain conditions are given by Miller 
(Appendix D of Reference 111). The rosulting material properties matrix, 
[H] , is 




l+v 


k-y,.. 

ii-2v 
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L symmetric 


I ( V ®11‘’22\ / 811«12\ / V »ll^33 
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(A-32) 


where B “ I ' (l+'O 

and tlie elemfsnts of the incremental dllatational strain vector are 




where Cq '= «AT + PAM 


r . > 


b.j 


(A-33) 


(A-3A) 
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AT and AM representing the Incremental changes in temperature and moisture 
content, respectively. Thus, for each matrix material element that has 
exceeded its elastic limit at a particular temperature and moisture con- 
dition, Equation (A-32) Is used to define Its material properties so that 
Its elemental stiffness matrix. Equation (A-9), can be formed and Incorpor- 
ated into the global stiffness matrix, [K], for the current load increment. 

An additional feature of the present analysis is that it is not re- 
stricted to monotonic loading. That is, unloading from the plastic region 
is possible for those elements that undergo stress relief due to local de- 
formation or crack formation. It is assumed that the material unloads 
linearly elastically, with a modulus equal to the elastic modulus. If 
total unloading were to occur, the total plastic strain sustained by the 
element up to the point of unloading would be retained. If reloading sub- 
sequently occurs, the material follows the stress-strain response described 
by the elastic modulus until it attains the stress level from which It 
began to unload. For any loading past this point, the material continues 
to move along the original elastoplastlc curve. 

A. 5. 2 HYGROTHERMAL DEPENDENCE OF MATERIAL PROPERTIES 

As noted previously, the alteration of the properties of the consti- 
tuent materials of a composite as they are exposed to a changing environ- 
ment has a significant effect on the microstress state of the composite. 

The present analysis scheme evaluates the stress-strain behavior, the 
coefficient of thermal expansion, the coefficient of moisture dilatation, 
the yield strength, the ultimate strength, and the ultimate strain for 
each element at the beginning of each load Increment. 

Many past elastoplastlc analyses have input the entire stress-strain 
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response of the material under study as data tables, utilizing considerable 
amounts of core storage or peripheral storage in the process. If one were 
to do this for each stress-strain curve corresponding to a range of temp- 
eratures and moisture levels, a prohibitive amount of storage would be 
required, and even on today's larger computers, this would significantly 
limit the size of model that could be analyzed. The present analysis uses 
a three-parameter model to parametrically describe the stress-strain curve 
of a material at any temperat'ire and moisture of Interest. The model used 


ij based on tlie Rlchard-BlacU lock equation [32], vWilch is capable of 
acojrately modeling material behavior which exhibits large amounts of 
plastic deformation. The general form of the Richard-Blacklock equation 



(A-35) 


where E is the elastic modulus of the material and Op and n are independ- 
ent parameters as Illustrated in Figure A. 2. 



FIGURE A. 2 

Richard-Blacklock Representation of Stress-Strain Curves 
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For the present analysis, Equation (2.35) is modified to describe an 
octahedral shear stress-octahedral plastic shear strain curve, as the 
primary material failure criterion employed is the ultimate octahedral 
shear stress criterion. Fortunately, uniaxial test data for a given 
material can be converted directly to octahedral shear stress and strain. 
The octahedral form of Equation (A. 35) is 



(A-36) 


where to is octahedral shear stress as defined by Equation (A-28), and 
To end n are again Independent parameters which fit the curve to empirical 
data by means of a numerical least squares curve fit procedure (p. 76 of 
Reference [1]). The term e is the octahedral shear strain, defined as 

7 - (1/3 (A-37) 

and E Is the slope of the initial linear portion of the octahedral shear 
stress-octahedral shear strain curve, being related to the elastic modulus, 

E, as follows 

® w-38) 

The tangent modulus of the octahedral shear stress-plastic octahedral 
shear strain curve, 2M'p, is the quantity required in forming the material 
properties matrix, [H], shown in Equation (A-32). The tangent modulus, 

2Mx, can be related to the tangent modulus, E-p, of the octahedral shear stress- 
octahedral shear strain curve by 

2Mt = -P5- (A-39) 

E-Et 
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where Ef can be obtained by differentiating Equation (A-36) with reapect 
to e. This yields 


K- I 
Ee 

1 + -- 
- 5o 


1/n 


(A-40) 


Thus, for a given octahedral shear strain, the tangent modulus, 2Mj., 
can readily be calculated. It has been found that small changes of tem- 
perature or moisture do not drastically alter the stress-strain curve of 
a matrix material, but rather modify the material properties in a uniform 
manner. In other words, a functional relationship between the parameters 
Og, n and the elastic modulus, E, and temperature and moisture states can 
be found. This is done by fitting Equation (A-36), In a least squares 
sense, to octahedral shear stress-octahedral shear strain curves obtained 
from tensile test results in various environments. A somewhat more de- 
tailed description of this procedure is available on pp. 73-78 of 
Reference [IJ. 


